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ABSTRACT 

Many of the recently developed high-resolution schemes for hyperbolic conservation 
laws are based on upwind differencing. The building block of these schemes is the averaging 
of an approximate Godunov solver; its time consuming part involves the field-by-field 
decomposition which is required in order to identify the “direction of the wind.” Instead, 
we propose to use as a building block the more robust Lax-Friedrichs (LxF) solver. The 
main advantage is simplicity: no Riemann problems are solved and hence field-by-field 
decompositions are avoided. The main disadvantage is the excessive numerical viscosity 
typical to the LxF solver. We compensate for it by using high-resolution MUSCL-type 
interpolants. Numerical experiments show tha t the quality of the results obtained by such 
convenient central differencing is comparable with those of the upwind schemes. 
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INTRODUCTION 


In this paper we present a family of non-oscillatory, second order, central difference 
approximations to non-linear systems of hyperbolic conservation laws. These approxima- 
tions can be viewed as natural extensions of the first-order Lax-Friedrichs (LxF) scheme 
In particular, total-variation and entropy estimates are provided in the scalar case, an 
unlike the upwind framework, no Riemann problems need to be solved in the case of sys- 
tems of conservation laws. The use of second-order piecewise-linear approximants instead 
of the first-order piecewise-constant ones, compensates for the excessive LxF viscosity, and 
results in second-order resolution Riemann-soh er-free family of central difference schemes. 

The paper is organized as follows. In Section 2, we derive our family of high resolution 
central differencing schemes, using the LxF solver together with MUSCL-type interpolants. 
Thus, at each time-level we reconstruct from the piecewise constant numerical data, a non- 
oscillatory piecewise linear approximation of second order accuracy. We then follow the 
evolving solution to the next time level, and end up by projecting it back to a piecewise 
constant solution. The result is a family of schemes which takes an easily implemented 
predictor-corrector form. The resolution of our method hinges upon the choice of certain 
local numerical derivatives with which one reconstructs the piecewise-linear MUSCL-type 

interpolants from the piecewise-constant data 

In Section 3, we concentrate on the scala - conservation law. We discuss a variety of 
choices for numerical derivatives, and prove that the resulting scalar family of schemes, un- 
der the appropriate CFL limitation, satisfies both the Total Variation Diminishing (TVD) 
property and a cell entropy inequality. These properties guarantee the convergence to the 
unique entropy solution, at least in the genuinely non-linear scalar case. 

In Section 4, we describe several ways to extend our scalar family of central differ- 
encing schemes to systems of conservation laws. The main issue lies again in the choice 
of vectors of numerical derivatives. First, we describe a component-wise extension for 
the definition of these vectors, which share the simplicity of the scalar family of schemes. 
Next, we demonstrate the flexibility of our central differencing framework, which enables 
us to incorporate characteristic information-whenever available, into the definition of nu- 
merical derivatives. We continue, by using :his characteristic-wise framework to isolate 
the contact wave where the Artificial Compression Method (ACM) is employed, while 
treating the more robust sound waves using :he less expensive component-wise approach. 
We end up by presenting a corrective type ACM, which is implemented in a component- 
wise manner. This both improves the contact resolution, and retains the simplicity of the 

Riemann-solver-free scalar approach. 

Finally, in Section 5 we present numerical experiments with our high-resolution non- 


L 


oscillatory central difference schemes, and compare the results with the corresponding 
upwind-based ones. 

Both the quantitative and qualitative results for a representative sample of compress- 
ible flow problems governed by the Euler equations, are found to be in complete agreement 
with the resolution expected by the scalar analysis. Taking into account the ease of im- 
plementation, robustness and time performance, these results compare favorably with the 
results obtained by the corresponding upwind-based schemes. 
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A rAMILY OF HIGH-RESOLUTION CENTRAL DIFFERENCING METHODS 

Many of the recently developed high-resolution schemes, which approximate the one 
dimensional system of conservation laws 

Tt + = °’ <2 ' 

are based on upwind differencing. The prototype of such upwind approximationsisthe 
Godunov sCrf it computes a piecewise constant approximate solut.on over ceils 

width Ax = x, + i - ly-J. which is of the form ’ 

mx,t)=Vj{t), <*<x i+ .. (2 ' 2) 

To proceed in time, the Godunov scheme first evolves the piecewise constant , sohrUon, 
w(x t) for a sufficiently small time step At. Initiated with e(x,t), equal ( • ) 

£ successive sequence of non-interacting Riemaun problems. Their resuUing soiut.on at 
time level t + At, can be expressed in terms of the R.emann solver, R[j, 

u ( x,t + a.) = .iW. ‘hW). ** * 1 s (2 ' 3) 

This solution is then projected back into the space of piecewise constant gridfunctions, see 
Fig. 2-1, 


„( t + A.) .«-(.,* + At j^/^vly.tEAOdy, 


+ AO 



Fig. 2-1 
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Integration of (2.1) over a typical cell [*,_j,s, +} ] x |t,t + At] yields 

■V(t + At) = u,(t) + A|/(*(0+; vy-.ft), „,.(!)) - /(*((,+; „ iW) „. +l(())Jj A s At (2 

IS shows the upwind property of the Godunov scheme Namely if the h A ** ■ • 
speeds throughout the relevant neighbouring cells I* * , , 

tively, negative), then (2.5) is simplified into „,« +' <TZ ,Z 
(respectively, Vj(t + At) = vAt) — Af f[v (t\\ t( wm ' MMO) /(»/-i(0)] 

«... J. f 

2 r - - - “ 

in^Zr ( ° r aPPrOXimate) S0,Uti0n ° f the Z”“l e , R alnn n d 

intenlrsCZrr*" flC ' d ' by ' Wd Proposed by Roe [19], which 

Instead, in this section we propose a high resolution approximation of f 2 1 1 v n • 

based on the staggered form of the Lax-Friedrichs (LxF) scheme, ' ’ “ 

”>>j(‘ + A ‘) = j[n ) ' + uy„]-A[/(u ) . +1 (t))_ /( „. (())1 (26) 

a g^nrsimpiicityZ'er'the'Zwhn^GZinov whe^e^^sJ^We^^erv^Uw.I^'c^^ 


oo 

v j+i(t + At) = v(x, f + At) = -L r i+l R t 

Ax Jx . K 


x ~ x j+i 

At ; V i» U J+l)dx, Xj < x < x y+1 . (2.7) 



+ 


f 


*» ^l+J *1+1 


Fig. 2-2 
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U V * f thP T xF scheme (2.7) , stems from the fact that unlike the Godunov 

The robustness of the Lxfc scneme, ^ ), ,1 left- 

is evident from the viscous form [23] 

Vj(t + At) = v,W - iM/(« J+ i(0 - /(®<->W)l+ (2.8) 

+M<2, + i A W‘> - (01. A W') 2 ”> + ‘ W “ ” ,W ' 

Indeed, the class of upwind schemes is characterised by a numerical viscosity ' 

• ■ n'j l| A il {here A + . refers to an approximate average of the Jacobian o 

c^fcL.l X |M+Atl, e,„ >22,). By the CFL .imitation, this amount 

- numericai viscosity is always iess than the amount of numerical v.scos. y presen Un^ 

central LxF scheme, whose non-staggered form corresponds to Q - ■ 

the upwind Godunov-like approximations have better resolution than the central Lx 

approximation, though they both beiong to the same Cass of 

This is one of the main motivations for using upwind schemes as bui 

modern shock capturing methods of higher (than first-order) resolution, e.g. [ ). I U 

Alternatively, our proposed method will use the simpler central LxF solver as the 
building block for a family of high-resolution schemes. In this manner we shall retain th 
LxF main advantage of simplicity: no Riemann problems are solved and hence ft - y 
field decompositions are avoided. The main disadvantage of excessive numerical v.scos. y 
will be compensated by using high-resolution MUSCL interpolants, (24], instead 
first-order piecewise constant ones in (2.2). 

"his end, at each time level we first reconstruct from (2.2) a piecewise linear ap- 
proximation of the form 

z*«. •)-.,(•) + (*-*>£•}• *i- t **<*,♦*• (29a) 

This form retains conservation, i.e., (here the overbar denotes the [*y.*f+il ’ cel1 average), 

Lj{x,t) = v{x,t) = Vj{t)\ 

second-order accuracy is guaranteed if the so-called vector of numerical derivative, £* 
which is yet to be determined, satisfies 

d 


J_v' 

Ai 3 


dx 


■v(x = Xj,t) + 0( Ax). 


(2.96) 
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Next, we continue with a second stage, similar to the construction of the central LxF 
recipe: we evolve the piecewise linear interpolant, (2.9), which is governed by the solution 
o successive sequences of noninteracting Generalized Riemann (GR) problems, (1], see 


v(x,t + At) = GR(x,t + At; Lj{x, t), L j+1 (x, <)), 


X; < X < x 


i+1* 



X 


Fig. 2-3 

Finally, the resulting solution is projected back into the space of staggered piecewise 
constant gndfunctions 

u, + . (t + At) = tr(*, t + At) ZE j- j‘‘*' „(„, , + Af)dji, z ( < lS XH1 
In view of the conservation law (2.1), the last integral equals 

f x . 7 L i ( x ’ ^ dx + f x L j+1 ( x , t)d. 


( 2 . 10 ) 


»,■+*(* + Af) = J_ 
2 Ax 




1 f t+At ft+At ] 

A x[Jr=t f( V ( X i+l> T )) dT - J T=t f{v{Xj,T))dT\. 


( 2 . 11 ) 


The first two linear integrands on the right of (2.11) L (x t) and T (*• t\ 

a u • x 1 J ’ and Lj +1 [x,t), are given 

by (2.9a) and can be integrated exactly. Moreover, if the CFL condition 




( 2 . 12 ) 
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is met then the last two integrands on the right of (2.11), /(v(z, ,r}) and f(v( ,+i> ))> 
are smooth fnnctions of r; hence they can be integrated approximately by the midpoin 
rule at the expense of O(At) 3 local truncation error. Thus we arrive at 

u,. + .(t + At) = i h (.) + », •«(»)) + + T» - /W *'* + l”', 


By Taylor expansion and the conservation law (-.1), 


t>(iy, t + = Vj{t) - ^A/', 


(2.14) 


. *i i + £1.} within the permissible second order 

mav serve as our approximate midvalue, t 2 J, . r , 

may serve as PP approximate numerical derivative of the 

accuracy requirement. Here, stanas ior ai. 

flux f(v(x = Xj,t)), 

-/'. = — f{v{x = a:y,t)) + O(Az), 


(2.15) 


Ax'* 

which is yet to be specified. 

We should emphasize that while using the central type LxF solver, we integrated over 
the entire Riemann fan, see .(** + At) in (2.10), which consists of both the left- and 
rightgoing waves. On the one hand, this enabled us to ignore any detailed k n° wle ^ e 
the exact (or approximate) generalized Riemann solver GR(-; -, -); on t e ot er an 
enables us to accurately compute the numerical flux, £*“/<•<*•'»*. wh~ values are 
extracted from the smooth interface of two non-interacting Riemann problems. 

In summary, our family of central differencing schemes takes the easily implemented 

predictor-corrector form, ^ + ^ , (() _ 1 A/ , (2.16a) 

» j+i (i+A«) = IhW +^.W)+ “ ' ( * ( ‘ + T H1 - (2 ' 16!,) 

Here the numerical derivatives of both gridfunctions, {v,} and {/,}, should obey the ac- 
Imc constraints (2.9b) and (2.15). In this manner the second-order accurate corrector 
step (2.16b), augments the first-order accurate predictor step (2.16a), and results in a high- 
resolution second-order central difference approximation of (2.1). 


Remarks: 


1. The choice jjwJ- 


_ _wi = o in (2.16), recovers the original first-order accurate LxF 
Ai J i ' 


scheme (2.6). 


7 



2. If instead of (2.6) we use the non-staggered version of the LxF scheme, 

v } {t + At) = i(t , j+1 (t) + v, •_!(*)] - j[/(v, +1 (i)) - /(t ,,_,(*))], (2.17) 

and repeat the reconstruction, evolution and projection steps described above, then 
the resulting high resolution central differencing approximation amounts to 




*v(< + T ) = v,(t) - 1 a/', (2.18a) 


To guarantee the desired nonoscillatory property of these approximations, the two free 
ingredients at our disposal - the numerical derivatives and £/j, should be carefully 
chosen. This issue will be discussed in the next two sections. 


3. THE SCALAR PROBLEM 

In this section, we are concerned with non-osciliatory high- resolution central differenc- 
ing approximations of the scalar conservation law 

d u d 

~dt + ~di = °* (3.1) 

Our family of high-resolution central differencing schemes (2.16) can be rewritten in the 
form 

”'♦*(* + = jM<) + ®i«Wl - A[ft +1 - j,], (3.2a) 

where the so-called modified numerical flux, g jy [18], is given by 

9i = fM‘ + ~)) + in', „,(( + = „,(() - i A /;. (3.26) 

^ ere ’ a x v j an approximate slope at the grid point Xj, 

1 > d , 

Ax Vj dx V ( X ~ X >' + °( Ax )> (3.3a) 

and 'hfj 1S the numerical derivative of the gridfunction {/,}, 

1 d 

Ax^ = dx^ v ^ x = x i' 0 ) + O(Ax). (3.36) 

The constraints (3.3) with smooth (= Lipschitz continuous) first order perturbations on 
their right, guarantee the second-order accuracy of the central differencing schemes (3.2). 
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In order to ensure that these schemes are also non-oscillatory in the sense to be described 
below, our numerical derivatives, should satisfy for every gridfunction to - 

0 < w ) • sgn{ Au i± i) < Const. • \MinMod{Aw j+ i, Au^}!- ( 3 ' 4 °) 

Here, the MinMod{; •} stands for the usual limiter, 

MM{ i, y} = MinMod{x, y} = ^syn(x) + sgn{y)\ ■ Min{\x\ , |y|) , (3.46) 

and can be similarly extended to include more (than two) variables. The constraint (3.4) 
is required in order to guarantee the Total Variation Diminishing (TVD) property for the 
family of central differencing schemes (3.2). We recall that TVD is a desirable property 
in the current setup, for it implies no spurious oscillations in our approx.mate solut.on 

v(x,f), [7]. 

However, it is well known, e.g. [7], [18], that one cannot satisfy both the accuracy 
requirement, (3.3), and the TVD requirement, (3.4), at the non-sonic critical gndvalues, 
«■, where Av + . • Au,_i < 0 # a(«y). Therefore, the second-order accuracy requirement, 
(3.3), must be 'given up at these critical gridvalues. Difference schemes with (formal) 
second order of accuracy at all but these critical gridvalues may be classified as having 
second order resolution in the sense that the local truncation error is almost everywhere 
O(Az) 3 , and the overall second-order accuracy does not seem to be degraded in such cases, 

at least in the L^norm. 

We shall verify the TVD property of the central differencing schemes, (3.2), with the help 
of 

Lemma 3.1: The scheme (S.Sa) is TVD, if its modified numerical flux, g j} satisfies 
the following generalized CFL condition , 


Ao,.i 1 

Ay i+ f = g j+ i - Qi- 

Av i+| 2 


(3.5) 


Indeed, by (3.2a), the difference v j+ i{t + At) - «;-i(* + At) equals 


v i+ >(t + At) At) - Av i+ »(;,- A AUj . + P + AUi ^2 + A Av j _J‘ 


Condition (3.5) tells us that the terms inside the parenthesis are positive and TVD follows 
along the lines of [7], 

TV{v{t + At)) = ^2 K+i (* + At ) “ v i-^ + - TV ( V W' ^ 
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Equipped with lemma 3.1 we turn to 


Theorem 3.2: Let the numerical derivatives ^v'. and A-ff. in (3.3) be chosen such 
that the TVD requirement (3.4) holds, say, 

0 < u' ■ sgn{Av j±1 _ ) < Const „ • \MM{Av j+ 1 , Av y _i}|, Const v = a, (3.7a) 

0 < /; • sgn{Av j± i) < Const f ■ |MM{Av y+ ±, A Vj _ , }|. ( 3 . 7 &) 

Assume that the following CFL condition is satisfied 

A • max |a(v ; )| < /?, /? = \ ___£ < — -(y^ + 4a — a 2 — 2). (3.8) 

Then the family of high-resolution central differencing schemes (3.2), (3. 3) is TVD. 
Proof: By (3.2b) we have: 


A W~ l + 8' , A ^T 1 - 

< jlII ~ fM t + t )) | ,®/+i(* + ?At) - Vj {t + ^) , l , Av> } +i. 

— I i A t\ /. a # \ * " -4- — 1 3 


(3.9) 


;,(« + *) 


j+k 


8 At/.-. i '■ 

J+5 


Our CFL condition (3.8) implies that the first term on the right of (3.9) does not exceed 


, i +f))-/(v,-(t+ 4*)) . ^ . 
<'/-«(*+¥) -v,{t+ f) l ~ P ‘ 


(3.10) 


Using the midvalue vfft + f ) in (3.2b), we can estimate the second term on the right of 
(3.9), 


+ ¥)-«,•(* + ¥) . j + a , A/; + f 


Av 




2 'A V , 


I- A/; +t s /; +l _ yj. 


where in view of (3.7b) and (3.8), 


A /' , 

2 1 < max 


/; 


/;• 


1 Av 


;'+ 1 


J+7 


Av 


>+T 




< Constf < —aj3. 

A 


(3.11a) 


(3.116) 


Finally, the TVD requirement, (3.7a), gives us an upper bound for the third term on the 
right of (3.9), 


A -; + i 


- max(l a^t 1 ’ 'a^ I) s “- 


(3.12) 
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Using (3 


.10), (3.11) and (3.12), we find that (3.9) boils down to the quadratic inequality 


0(1 + ;“0) + 5“ 5 b 


whose solution yields the OFL limitation (3.8). 

Remarks : 

1. The values a which permit a positive solution of (3.8), 0 > 0, are 0 < a < 4. 

2 The TVD constraints (3.7) with a = 0 yields «} E 1', = 0, which recovers the 
staggered LxF scheme (2.6) with the corresponding CFL condition 0 - 5 - 

TUT . • 4.- (o :<5 a sufficient but not necessary condition for the TVD 

3. The CFL restriction [3.5) is a sumcieni uut ^ 

property. In practice one may use higher values of 0, up to 0 < 5. 

4. A similar analysis carried out for the non-staggered form, (2.18), yields 

/3 < —(Vi' + 4a - a 2 - 2) 
a 

instead of (3.8). In practice one may use 0 < 1 in this case. 

We shall now discuss a few examples of numerical derivatives, which -tain both th 
second order resolution constraint, (3.3), and the TVD constraints, (3.7). As 
example for the numerical derivative, vj, we choose 

v*. = MM{ Au J+ i , Aw,. 1}. (3.13a) 

v * ctroua “discontinuity”, where the order of accuracy is less 
This choice may oversmear a strong aisconLinumjr , 

significant. A preferable second choice, which allows for a steeper slope near such discon- 
tinuities and yet retains higher accuracy in smooth regions, is given y 

of = “ v i- *)> “ A, i- * 3 ' 136 * 

The limiting parameter a can range between the values a = 1, which cor - s P° nd ® , ° 
the basic MinMod limiter in (3.13a), and up to a < 4, which is permitted by the 
condition (3.8). Similarly, the flux numerical derivative may be chosen as 

(3.14a) 


/; = MM{A/ i+ »,A/,_.}, 


which is a special case of 


f. = ;-(/>+! 


(3.146) 
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A simpler alternative for (3.14) is given by 


fj = 


(3.15) 


»here »} is already computed by (3.13). We observe that this choice saves half the com- 
putation time of the MinMod operation; yet, it requires the computation of the Jacobian 
A[ Vj ), when dealing with systems of conservation laws. 

The numerical derivative chosen in (3.13a), (3.14a) satisfies (3.7) with a = 1 which 
mrphes the TVD property under the CFL limitation (3.8) with 0 = l (v /7- 2) „ 0 32 

* ~ 1 denVat ‘ VC chosen in < 3 - 13b ). (3-14b) clearly satisfies (3.7) and conse- 
TVD Pr ° Perty ' fW €Very P-nrissible a, 0 < « < 4. We summarize the above 

Corollary 3 . 3 : Let the numerical derivative be chosen by 
ht the flux numerical derivative be chosen either by 

fj = a ( v j)v'j. 


(3.16a) 


(3.166) 


or 


/; = mm{a/. + i,a/._i}. 


(3.16c) 

2ZfL iZZl h ' Sh rCSOht<0n C ‘ ntr °' diff ‘ ,mcins schcm “ < S S >- (fM) is TVD under 


A • maxj\a(vj) \ < 0, fl = 1(^7 - 2) « 0.32. 

Similarity W6 have 


Corollary 3 . 4 : Let the numerical derivative be chosen by 


Vj = MM{2Av j+ ^, ~(vj +1 - v ; _i),2 At/ y _i} ; 
let the flux numerical derivative be chosen either by 

(3.17a) 

fj = a ( v })v'j, 
or 

(3.176) 

/' = l(/ )+1 - /,_,), 2A/ y _,}. 

(3.17c) 
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Then the family of high resolution central differencing schemes (S.2), (3.17) is TVD, under 
the CFL condition, 


X • maxj\a(vj)\ < /?, 


P : ~ r (V2 — l) ~ 0.21. 

It 


Remarks : 

1 We note that the CFL limitations in Corollaries 3.3 and 3.4 are not sharp. In the 
first case, (3.16), where a limiter parameter a = 1 was used, the reconstruction step 
is a TVD operation; replacing the exact TVD evolution operator by the midpoint 
rule in (2.11) together with the final averaging step is also TVD, under the CFL 
limitation (3 < k- Similarly, one can argue that in the second case, (3.17), where 
a limiter parameter a = 2 was used, the averaging step retains the TVD property 
(though not necessarily the entropy condition), as long as the CFL condition /? < j 
is met. Indeed, this CFL condition was verified as the stabilitiy limitation, by t e 
numerical experiments reported in Section 5. 

2. Recently, non-oscillatory schemes were constructed, such that by sacrificing the TVD 
property, they achieve higher (than second-order) resolution including the critica 
gridvalues, e.g., the UNO scheme in [12] and the ENO class of approximations in 
[91 To implement such ideas within our framework, one can borrow their definition 
of numerical derivative. For example, instead of the TVD choices (3.4), our central 
differencing scheme (3.2) may be augmented by the UNO choice (here A vy - «y+i 

2 Vj + Vj—i), 

= MM{ + \MM{ AVi, AS), A» i+ * - A a vy, A*vy+i)}. (3.18) 

Theorem 3.2 and its corollaries 3.3 and 3.4 demonstrate high- resolution central differ- 
encing methods which satisfy the non- oscillatory TVD property, and hence are convergent 
to a limit solution «(*,«)• To guarantee that this limit solution is the unique entropy so- 
lution of the scalar conservation law (3.1), we shall appeal to the following cell entropy 

inequality, see [10], 

U{v j+ ^{t + At)) < ^[U{vj) + U(u J+ i)] - A[G i+ i - Gy]. (3T9) 

Here U{u) is a convex entropy function and G, = G(tfy+i, vy,vy-i) is the numerical entropy 
flux which is consistent with the corresponding differential one 


G(u,u,u) = F(u), F (u) = f («). 
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We recall that Lax has verified such cell entropy inequality for the LxF scheme, [14]. 
Following Lax, we will continuously deform Vj into v/+j. 


v(s) _ svj + (1 - a)t>y +1 , i>(0) = v y+1 ,v(l) = Vj , 

and m a similar manner, we will further deform v(s) into u y+1 , 

w(r,s) = ru( 5 ) + (1 - r)v j+1 , i>(0,s) = v i+u v(l,a) = «(«). 

In the Appendix we prove 


(3.20a) 


(3.206) 


Lemma 3.5: Let g(v) be a piecewise differentiable interpolant of the gridfunction {g } 
Then the following identity holds, 

U{v j+ i.{t + AO) = -\U{v j+l ) + U{ V j)} U'{v)g'(v)dv - R” +l _{g{v)). (3.21) 

Here the residual term, R j+i _(g) = i?^,( 5 ( w )), is given by> 

£ ffv"(v(r,s)) ■ [1 - Aj'(„(r, 5 ))] . [i + \g'(v{s))) is dr. (3.22) 

Adding and subtracting 

f v ’ + u '(u)f'(u)du = F( V j +1 ) - F( V j), 

then after integration by parts, the right hand side of (3.21) will amount to: 

U ( v 3+^( t + A *)) - + u ( v j)\ - A[F’(u y+1 ) - F’(vy)] - A U'(v) ■ ( ff (u) - /( v ))|JJy+i 

+X L’ +l • (s( v ) ~ f(v))dv - R^i(g(v)). 

Consequently, the inequality 

A f^' U"(v) ■ (g(p) - /(,))* - (,,(„)) < o, (3.23) 

provides us with a sufficient condition for the family of central differencing schemes (3 2) 
to satisfy the cell entropy inequality, (3.19), with numerical entropy Hurt G, = Flv) + 
U'M ■ ()M - /(»,)). This brings us to ’ 1 ’’ + 

Lemma 3.6: Let g(v) be the piecewise linear interpolant of the modified flax gridfunc- 

Uo n {gj}, 




Ay v i) + min(v y ,vy +1 ) < v < max(i/y, vy +1 ). 


(3.24) 
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Assume that the central differencing schemes (S.S), satisfy the TVD constraint (consult 

(S-7», 

0<u' r sgn(Av j±i ) < Const, -\MM{I\v^,,txv+_,}\, Const, = a < 1, (3.25o) 

Wh ‘ rC Au; + ,=Au )+r (l-A-(max/>W)'Au, + i) + r- ( 3 - 25t ) 

The entropy dissipative Liter in (S.SSb), is intrcduced in order to prevent the non expansive 

entropy violating rarefactions, consult (18, Section 8f. 

Moreover, assume that the flux numerical derivative satisfies the TVD constrain . 


Const f 


0 < /j • sgn(Av J± i) < Const f • \MM{Av i+ i, A • - P 


(3.25c) 


so that the CFL condition (8.8) holds. 
Then the following inequality holds 


A r‘(g(v) -/(»))* -i? t5 (sW)<0. U(u) = \u\ (3.26) 

Jv,- 


Remarks : 

1. We observe that in the Genuinely Non-L.near (GNL) case, where, say, /" > 0, the 
entropy entropy (3.25b) becomes effective only in rarefaction cells where At) J+ i > 0, 
in agreement with [18]. It retains the second-order resolution of the central differ- 
encing schemes (3.2) , except for a finite number of critical cells which contain strong 
rarefactions, (A V;+ ,) + ~ 1, where it reduces (3.2), (3.3) to the original LxF scheme. 

2. Lemma 3.6 applies to choices of numerical derivatives, subject to the TVD con- 
straint (3.7a) with 0 < a < 1. In practice, higher values, a > 1, can be used. 

Lemma 3.6 - which is proved in the appendix, shows that our central differencing TVD 
schemes (3.2), (3.7) fulfill the sufficient condition (3.23) and consequently the cell entropy 
inequality (3.19), with respect to the quadrat: c entropy function U(u) 2 u . 
limit solution of our central TVD schemes, u(:c,t), satisfies 

£(!«*) + -|-(F(tt)) < 0, F(u) = J uf'{u)du. 
dt K 2 dx J 

This singles out «(„,*) as the unique entropy solution of (3.1). at least in the GNL case 
[2]. We have shown 

Theorem 3.T: Consider the GNL scalar conservation law (S.l). It is approximated 
by the family of high resolution central differencing schemes (S.S), ( S S ) ^ h saUs/y 1 e 
TVD and entropy constraints, (S.SS). Then, if the CFL condition (S.S) holds, we have. 
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1. Second-order resolution; 

2. Total Variation Diminishing property; 

S. A consistent quadratic cell-entropy inequality; 

and, as a consequence of 2. and 3.: 

*■ thc >'"» "ntra! differencing scheme , converge to the unique physically r el- 

evant solution of the GNL conservation law (S.l). 


We shall close this section with some scalar numerical examples. We consider the 
approximate solution of the inviscid Burger’s equation 

u t + (~u 2 ) * = 0. (3.27) 

using several of the previously mentioned central differencing schemes. They include: 

1. The first-order LxF scheme in its non-staggered form (2.17). 

2. The second-order non-oscillatory central differencing scheme (2.18), (3.13a), (3.15). 

This is the ordinary non-staggered version of our central differencing which will be 
referred to as ORD. 

3. The second-order non-oscillatory central differencing scheme (3.2), (3.13a), (3.15). 
This is the st aggered version of our central differencing which will be referred to as 


Equation (3.27) is solved with two sets of initial conditions. In the first case, we have the 
smooth 1-periodic initial data, 

u(x,0) =m(n)i (3.28a) 

in the second case we consider the Rieman initial data: 


u(x, 0) = 


1 x < 0 
0 x > 0 


(3.286) 


The well known solution of (3.27), (3.28a), e.g. (15), develops a shock discontinuity at 
t c » 0.31. Table 3.1 shows us the L, norm of the errors at the pre-shock time t = 0 15 

It indicates the first order accuracy of the LxF scheme in contrast to the second order 
accuracy of our central differencing, ORD and STG. 

Next, in Table 3.2 we recorded the same L x errors at the post-shock time t = 0.4. The 
presence of a shock discontinuity reduces the global Z,, error to first order. However, the 
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central differencing STG scheme performs somewhat better than the central differencing 
ORD scheme and they both have better resolution than the first-order LxF scheme in 

shock-free zones. 

This behavior is amplified in the second case of solving the Rieman problem (3.27), 
(3.28b), see table 3.3 and Fig. 3.1. Once more, we observe that the STG scheme has 
somewhat better resolution then its non-staggered counterpart ORD. Yet, the CFL lim- 
itation in the non-staggered form, ft < 1, results in a better time preformance than the 
STG scheme which is subject to the CFL limitation ft < §. In either case, these easily 
implemented non-oscillatory central differencing outperform the first-order LxF one. 


17 



Table 3,1 : Li norm of the 


errors, Burger Equation, u(x,0) = sin( 7 rx), t = 0.15. 


Rel. 

STG 

Rel. 

ORD 

Rel. 

LxF 

NX 

3.70 

.000859 

3.92 

.002620 

1.93 

.023702 

40 

3.80 

.000232 

3.94 

.000667 

1.96 

.012249 

80 

3.81 

.000061 

3.93 

.000169 

1.97 

.006246 

160 


.000016 


.000043 


.003158 

320 


T able 3.2: Li norm of the errors, Burger Equation, u(x,0) = sin(Trx), t = 0.4. 


Rel. 

STG 

Rel. 

ORD 

Rel. 

LxF 

NX 

3.06 

.000849 

2.79 

.003612 

1.89 

.044449 

40 

2.82 

.000277 

2.59 

.001291 

2.06 

.023486 

80 

2.57 

.000098 

2.38 

.000498 

2.17 

.011383 

160 


.000038 


.000209 


.005235 

320 


Table 3.3: Li norm of the errors, Burger Equation, u(x,0) given in (3.28b), t = 0.3. 


Rel. 

STG 

Rel. 

ORD 

Rel. 

LxF 

NX 


.010921 


.014044 


.021956 

50 

1.99 

.005461 

1.99 

.007027 

1.94 

.011315 

100 

2.00 

.002730 

2.00 

.003513 

1.99 

.005685 

200 

2.14 

.001274 

2.00 

.001756 

1.99 

.002843 

400 
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4. SYSTEMS OF CONSERVATION LAWS 

In this section, we describe how to extend our scalar family of central differencing 
schemes to the one-dimensional system of conservation laws, 

du d 

at + = °- (4-i) 

Here u(x, t) is the unknown TV-vector of the form 

u = ( Ul (x,t),u 2 (x, t),...,ujv(x,t)) r , 
and /(u) is the flux vector, 


f ( U ) ~ (/i( u )»/2(u),...,/tv(u)) T , 

with an N x N Jacobian matrix, 



Our approximate solution at the gridpoint x, is given by the IV-vector 


v i ~ i v i, i) v i,2t • • • > v j,n) T i 

and the corresponding vector of differences, Av j+ > = Vj+1 - Vj , consists of V-components 
denoted by Av j+ i k = Vj+lik - v jtk . 

Equipped with these notations, our family of high-resolution central differencing schemes 
(3.2), (3.3), takes the form, 


V j+ i(t + At) = -\ Vj (t) + u y+1 (t)] - - ffi ], (4.2a) 

where the modified numerical flux, g } -, is given by 

9j = fMt + ^)) + -^v', Vj (t + y) = «/(*) ~ /;. (4.26) 

As before, the computation of Qj and Vj (t + f ) requires the numerical derivatives of the 

gridfunctions {v^} and {fj}. This time we have to choose two N-vectors of numerical 
derivatives, 


Ax 


v 'j = 


(4.3a) 


w J> 1,j i'.2» ’Jj.N) • (4.36) 

In the rest of this section, we shall describe the pros and cons of several choices for these 
vectors of numerical derivatives. 
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(4.4c) 


Our first choice is a component-wise extension of the scalar definition m Section 

this end we may use either (4.4a), . 

i 1.-1 N (4.4a) 

v ;. k = MM{Av i+ i it ,Av^}- k-i,..., . 

or the more general (4.4b), 

1, „ \ a & v . , k = 1,...,N, ( 4 - 4b ) 

t,' fc = MM{« jfa+i,* " ,_5 ' ' 

or instead, use the UNO-like numerical derivative in (3.18), 

= MM( Au j .,, 1 +|MM(A , v j -i.iiA 2 v j , k ), (4.4,.) 

A» )+ l, t - |MA/(A 3 U;,k, AS+u)}, * = x > ' ' • ’ 

A possible choice for the vector of numerical Hu e derivative may be ^ 

» — — 

This multiplication may be avo.de For example, we may use 

of numerical flux derivatives, in analog, to 13.14). For examp 

/j k = MM{A/ J+ ^n 

or alternatively, 




( 4 . 66 ) 


v. m f 1F E" (4 4 ) (4.6) avoids the use of the Jacobian 
the - *• — > 

‘"“he resulting central differencing schemes, 

wise definition of the numerical deriva ^ and consequently characteristic 

framework. Namely, no Riemann pro betwee „ the left and rightgoing waves 

decompositions-required in or er o is central differencing approach 

• “ 

available, in order to achieve improved teaolution information into the 

Our next choice showshowtoincorpor^ the^chara ^ ^ ^ ^ ^ 

definition of the numerical deriva . _ A satisfying, e.g. Ill], [19), 

A t = A(v h v„ 1), namely, an averaged Jacobian, A )+ ,, 

/(u,+i) - /(",) = ^,+i ' f”'-*' “ 
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and having complete real eigensystem {a. > R . , \ u _ i A7 

vector of differences Av . , onto Iff li *’ 2 ’ N - Let us project the 

Vi onto {* y+ i}, i.e. we use the characteristic decomposition 


iVj+ * • ^i+ijt, A: = l,... ,JV j 


(4.8a) 


where 


(4.86) 


(4.9) 


S + i* i t • Au y+ ^ , L r Rj =zS ijf k = 1, . . . , iV. 

Then the corresponding projection of the hurt vector of differences is given by 

Now, a possible characteristic-wise choice fnr the „ . , , 

M ’ - b ' ^ ^ averaged ^ ^ 

and the numerical flux derivatives can be calculated as 


(4.10) 


(4.11) 


j! f 

of numerical derivatives, if inste^of | Vt “°' 

/;,r = EM(s /t)ll s, t |, 1 ,Virt- i A- 

As an example, let us consider the Euler equations, 


(4.12) 


<9 

’ /=> 

d 

m 

dt 

m 

> . 

dx 

pu 2 

. u{E + p) _ 


-o, p=( 1-1)-{E-I pu *y 


(4.13) 


Here 


-a, :r ~ 

.here the 


X i+Lk are given by 


a 7+i.l = U.u.1 - 

and the right eigenvectors 


J+j,i u y+ i c , i a i „ = 7/ , ^ 

2 ->+ a J+ 1, 2 - V , , a y+ , 3 = U y+ 1 + c, u 


(4.14) 




1 

u - c 

L H — uc 


are given by 


R } + ^,2 — 
J+k 


1 

U 

L W 


R i+h < 3 ~ 




1 

u + a 

ff + uc 


(4.15) 


i+i 
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The average quantities on the right of (4.14)-(4.15) given in [19] are, 

< v?u> k _ <ypH> . = r 1)( a _ I a ,), 

< ,/P > < > V 2 


Ej + P; 

Pi 


(4.16) 


where < w >= l(w, + t» J+1 ) denotes the used arithmetic mean. This brings us to the 
characteristic decomposition (4.8), where the characteristic projections, 

1 


&, + M = slot - rn), ^, = >1* r-” 1 ’ “<+i’ s= 2^‘ + r> ^’ 


(4.17a) 


(4.176) 

(4.17c) 


are expressed in terms of which are given by 

Hi = (7 - 1) • (E j+1 _ + Ih+S+k ' 

m = (™;+i - ^+i“i+|)/ £ i+£- 
We note that the second contact field associated with ^ + i, 2 » independent of the 
square root which is required only in the computation of the mean value sound speed 
c- Since this field is a linearly degenerate, it lacks the strong entropy enforcement 
typical to the other two genuinely non-linear field, and therefore, is usually smeared by 
numerical schemes. In our next choice of numerical derivatives, we incorporate only partial 
characteristic information. Namely, we isolate the less expensive (be., square root-free) 
characteristic projection on the contact field , and use the component-wise approach for 

the other two fields. 


Thus, we first separate the contact field, 




&Pi+k 


' 1 


= 

Am j+^ 

- • 

U 

. A E j+ i . 


. - 


W J 


(4.18) 


i+k 


and then define the vector of numerical derivative as 


r rl ‘ 


' 1 


Ap J+ £,Ap,_i 

Pj 

m' 

= MM{a j+ i 2 ,a _i} • 

< u > 

- MM 

Am^Am^i 

UJ 

J ' 2 » • 1 

A 

to 

V 

i 

. &E j+ i,&Ej_i 


(4.19) 


Similarly, computing the numerical flux derivative with a characteristic approach applied 
only to the isolated contact wave, 


' Af, + t. i ' 


A fj+k> 1 


■ 1 

A li + k. * 

= 

A fj+k’ 2 

-- "y+i,2 ■ a j+i,2 • 

u 

1 A 2 

. A fj+5< 3 - 


J+ j. 3 . 


2 U J 


(4.20) 


i+k 
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amounts to 


Si’ 1 1 1 A/y_l j 

/• 2 =MM(a ) . + , s a ) . + , 2) a J .., 2 a / .: 2 ). <«> -MM A?. + 1 2 ,A/. !’ 

J ^ *]> a f a / I 

2 J > A /j+i,3.A/ J _. >3 J 

(4.21) 

The latter approach enables us to use effectively the Artificial Compression Method 

(ACM) on the isolated contact field, e.g. [6], [7]. To this end, the contact wave isolated in 
(4.19) is modified by 


(4.22a) 


I 

= [MM{a J+ i 2 , dfy.i j} + OjVj) • 

“ 1 


Ap + 1, Ap_ 1 ] 

m 3 

XL 

1 *2 

-MM 

Am +i , Am - 1 

L J -1 


~u 

2 J 

i 

.A£,. + ,,A J 


where 0 } - and r 2 - are given by 


„ <* •+ i 2 - &_•_ i » I 

Q. ~ 3_ 2> 2 ' 


(4.226) 


a j+^,2i + i«i-» >2 r 

T} ~ MM {^ 1 ~ Xu i+^ 2 ’ a i+^,2> 2 ^ ~ 2} • (4.22c) 

Finally, we shall mention an alternative approach to the characteristic implementation 
of the ACM in (4.22). To this end, the Artificial Compression is implemented as a further 
corrector step to the component- wise approach presented in (4.2a) -(4.2b). This corrective 
type ACM takes the form, 

v-(t+ AO =*,(* + At) -«(W y+i 0< £ <1. (4.23a) 

Here, the compression coefficient, e, and W, are given by 
wj, A Wy+ , • Ai/,. + , > 0, 

u>j = MM{Avj_^,v RL ,Av j+ i}, (4.236) 

l w i+i> Aw i+$ • i < 0, 

where v RL is related to subcell resolution information (Harten, private communication 

[ 8 ]), 


W i + k = < 


VRL \v j+l (t+At) Vj^t+At) Ax y+ i(6 ;+1 +6 J _ 1 )| > = MM{Av } - _± ,Au y+ i} .( 4 .23c) 

The result is the central differencing scheme (4.2), appended by the component-wise 
definitions of numerical derivatives in (4.10) - (4.12), and complemented by the ACM 
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corrector step (4.23). This scheme, unlike the characteristic-wise implementation of the 
ACM in (4.22), enjoys the simplicity of the component-wise approach, and at the same 
time, enables us to deal effectively with the delicate contact wave. We remark that one 
should be careful not to overcompress discontinuities using such corrective type Artificial 
Compression: it should be implemented after th e rarefaction waves have evolved using an 
appropriately chosen compression coefficient e. 


5. NUMERICAL EXAMPLES 


In this section, we will present numerical examples which demonstrate the performance 
of our family of high resolution central differencing schemes for systems of conservation 
laws. We consider the approximate solution of the Euler equations of gasdynamics, see 

section 4, 


d_ 

dt 




m 

pu 2 

u(E + p ) 


= 0, P = il - 1) • {E 



m = pu. 


(5.1) 


We experiment with the following members from our family of high- resolution central 
differencing schemes: 


1. The central differencing scheme (4.2), (4.4a), (4.5). This is the component-wise 
extension of the scalar STG scheme presented in Section 3 and is therefore referred 
to by the same abbreviation. 

2. The central differencing scheme (4.2), (4.4b), (4.5) with a limiter value a = 2. This 
scheme is referred to as STG2. 

3. The component-wise UNO-type version of our scheme, (4.2), (4.4c), (4.5). It is 
referred to as STGU. 

4. The scheme STG with the addition of the corrective type ACM described by (4.23) 
is referred to as STGC. 

All the above examples use component-wise definitions for the vectors of numerical 
derivatives, and are based on the staggered grid formulation. Our last example is 
based on non-staggered LxF scheme, namely, 

5. The central differencing scheme (2.18), (4.4a), (4.5). This is the component-wise 
extension of the scalar ORD scheme pres ented in section 3 and is therefore referred 
to by the same abbreviation. 
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For the purpose of performance comparison we include here the results of several well 
known upwind and central schemes as well. These schemes include: 

1. The first order central non-staggered LxF scheme, (2.17), [13]. 

2. The first order accurate Godunov- type scheme of Roe, e.g. [7]. 

3. Harten’s second order accurate upwind ULTl scheme, [7]. 

4. Harten’s second order accurate upwind ULT1C scheme, [7], where Artificial Com- 
pression is added to ULTl in the linearily degenerate contact field. It is referred to 
as ULTC. 


We solve the system (5.1) with three sets of initial conditions. Our first example is the 
Riemann problem proposed by Sod [21] (abbreviated hereafter as RIMl), which consists 
of initial data 


„(*,<>) = \ vi ' * <0 ’ ^S 1 ’ 0 ’ 2 - 5 ) 1, f52) 

\ V r , x>0, IV = (0.125, 0,0.25) r . 

Table 5.1 shows the L\ norm of the errors. Though the results are field dependent, the 
“quantitative picture” is favourable with the central differencing schemes. Table 5.2 shows 
the time performance of the various schemes. All the schemes have time performances 
of order O(NX) 2 , where NX is the number of spatial cells. Figures 5. 1-5.4 include a 
comparison between the numerical solution and the exact solution (shown by the solid 
line), e.g. [3], [20], at t — 0.1644. As expected, the overall resolution of the first order 
schemes is outperformed by the second order schemes. 


We observe that our second order staggered schemes, STG, STG2, and STGU, and 
similarily, the second order upwind ULTl scheme, smear the shock discontinuity over two 
cells. The contact discontinuity, however, is more delicate: here we observe smearing of 
about 5-6 cells by the second order schemes, both in the central and upwind cases. We can 
also observe the over- and undershoots generated by both the upwind ULTl and central 
ORD. These unsatisfactory results suggest to introduce ACM in the contact field. For 
this purpose we present the upwind ULTC scheme and the central component-wise STGC 
scheme in Fig. 5.4. We note that the ACM is applied in STGC only at the last 10% of the 
time steps with the compression coefficient e = 0.625. This results in 2 cells resolution of 
the contact wave, and somewhat better resolution in the other waves as well. Yet, small 
over- and undershoots which are due to overcompression, still remain. 


Our second Riemann test problem (abbreviated hereafter as RIM2), is the one proposed 
by Lax [5]. It is initiated with, 

tV = (0.445, 0.311, 8.928) t , v r = (0.5,0, 1.4275) T , (5.3) 
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, f j • t^o- c c eg The density profile in RIM2 lacks 
and the results at t = 0.16 can be found m Fig. 5.5 5.S. me aen y y w 

the monotonicity wo had in RIM1, and therefore, it is more difficult for 

numerical schemes to recover the contact wave and the intermediate plateau winch 

follows. Consequently, the upwind schemes pe rform here somewhat better than the cent 

schemes: ULTC resolution is better than STGC which has more over- and ™dersh°o 

than before. We note that STG2 has better resolution and L, errors than STGU 

Selds. This is due to the fact that STG2 has steeper slope near discontinuities, consu 

S "y, the results of the nonstaggered central difference scheme ORD for both RIM1 
and RIM2 problems are presented in Figure 5.9. We recall that the CFL l.mitat.on in 
the staggered case, ff < ], is now doubled to be 0 < 1, consult Section 3 Moreover, a 
component-wise reconstruction of the vector of numerical derivatives, enabled us to avoid 
any Riemann solver in this nonstaggered case. Consequently, the ORD sc erne is wice 
faster than the staggered central versions based on STG, as well as the upwind scheme 
ULT1 which necessitates the (approximate) solution of a Riemann problem at each cell 
However, the resolution of this nonstaggered version, ORD, deteriorates, when compared 
to the staggered versions and the upwind methods. 


Our third problem, discussed by Woodward-Collela in [25], consists of initial-data, 


u(l,0) = I 


V t 0 < x < 0.1, 
v m 0.1 < x < 0.9, 
v r 0.9 < x < 1, 


(5.4) 


where « = „ = P. = 1, m, = m„ = m, = 0, * = 100, p„ = 0.01, p, = 100. A solid wall 
boundary conditions (reflection) is applied to both ends. The results are compared with 
the fourth order ENO scheme [9], in Fig. 5.10-5.12.’ The continous line is the result o 
the ENO scheme with 800 cells. We present the results of STG2 and ULT1 with 400 cells 
in Fig. 5.13-5.15 at t = 0.01, t = 0.03, and t = 0.038 respectively. We observe that the 
upgrade from the first order LxF scheme to the second order STG2, results in a >ub>tan .a 
improvement of resolution, see Fig. 5.10-5.15; moreover, STG2 compares favourably with 

the second order upwind ULT1 scheme. f cTr^wms 

In summary, we may conclude that when strong discontinuities are present, STG2 seems 
to offer the best results, STGC can be tuned to obtain sharp resolution at the expense 
of overcompression, and ORD version wee found to be the most economica . ur 
extensive numerical experiments done alorg these lines are reported in 116]. 

2 We t h an k A. Harten for allowing us to use his ENO results in [9]. 
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— ble 5,1 : Computation time of Riemann problems, results at t = 1.0. 


ULT1/C STG ROE 

RIM1 

1.23 
4.93 
19.81 


ORD LxF STGC 

1.43 
5.67 
22.74 


STGU STG2 NX 

1.47 
5.88 
23.49 


RIM2 

2.87 

11.54 

46.34 


1.23 0.74 

4.75 2.92 

19.32 11.68 


2.74 1.72 

10.93 6.83 

43.50 27.27 


0.69 0.22 

2.71 0.85 

10.71 3.37 


1.55 0.48 

6.16 1.90 
24.40 7.52 


3.24 3.35 

12.88 13.30 

51.46 53.20 


1.37 50 

5.43 100 

21.66 200 


3.07 50 

12.22 100 
48.83 200 


Remarks : 

1. Due to our method of implementation, ULT1 and ULTC have the same computation 
time. In fact ULTl is somewhat faster then ULTC. 

2. All the above schemes use a CFL number of 0.95, except for the versions, STG*, 
which use a CFL number of 0.475. 
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Table 5.2: Riemann problems, L i norm errors. 



Density 



Velocity 



Pressure 


Nx 

50 

100 

200 

50 

100 

200 

50 

100 

200 

Scheme 



Rieman Problem - RIM1, t — 0.1644 



LxF 

.03121 

.02460 

.01769 

.06651 

.04583 

.02814 

.03602 

.02458 

.01582 

ROE 

.01918 

.01308 

.00836 

.03224 

.02090 

.01145 

.01762 

.01109 

.00666 

ORD 

.01868 

.01026 

.00578 

.03315 

.01807 

.00959 

.01630 

.00861 

.00460 

STG 

.01495 

.00741 

.00409 

.02812 

.01105 

.00550 

.01232 

.00581 

.00294 

ULT1 

.01338 

.00806 

.00437 

.02933 

.01177 

.00820 

.01285 

.00736 

.00362 

STG2 

.01241 

.00619 

.00297 

.02449 

.01132 

.00494 

.01019 

.00487 

.00228 

STGU 

.01146 

.00544 

.00291 

.02300 

.00816 

.00403 

.00961 

.00432 

.00216 

STGC 

.00982 

.00322 

.00172 

.01994 

.00481 

.00276 

.00705 

.00270 

.00153 

ULTC 

.01269 

.00715 

.00361 

.02923 

.01761 

.00804 

.01283 

.00735 

.00362 




Rieman Problem • 

RIM2, t 

= 0.16 




LxF 

.12162 

.09044 

.06165 

.13523 

.09294 

.05557 

.15860 

.10767 

.06537 

ROE 

.06630 

.04334 

.02827 

.07397 

.04144 

.02192 

.08399 

.04826 

.02655 

ORD 

.06791 

.03824 

.02231 

.07158 

.03623 

.01709 

.07836 

.04056 

.01995 

STG 

.04972 

.02903 

.01776 

.04392 

.02416 

.01307 

.05118 

.02669 

.01426 

ULT1 

.04518 

.03572 

.01477 

.05570 

.02603 

.01094 

.06075 

.02841 

.01206 

STG2 

.03473 

.02129 

.01151 

.03369 

.01655 

.00849 

.03956 

.02037 

.00988 

STGU 

.03668 

.02152 

.01302 

.03323 

.01657 

.01046 

.03907 

.02031 

.01121 

STGC 

.02764 

.01291 

.00647 

.02285 

.01356 

.00836 

.02355 

.01409 

.00873 

ULTC 

.03001 

.01566 

.00872 

.05504 

.02545 

.01074 

.05997 

.02784 

.01183 


Remarks : 

1. All the above schemes use CFL number of 0.95, except for the staggered versions 
STG*, which use a CFL number of 0.175. 

2. The underlined results indicate the smallest Li norm errors in every column. 
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APPENDIX: ON A CELL ENTROPY INEQUALITY 

In this section, we provide the promised proofs for Lemmatta 3.7 and 3.8, whl ^ h Ver ' 

rt f nr our scalar family of high-resolution central difference 
ify the cell entropy inequality for our scalar ta y & 

methods. 

We begin with a proof of Lemma 3.7. 

Let R . i(fif) denote the difference, 

J ■ a 


B iH M - W.) + PWl - "WM* - u{v ^ {t + At))) - 


(A. 1) 


We now continuously deform ,(.) E .(.) = «, + P ~ +*»• “ " = ^ 

„. +I = „(1), see (3.20a). With this in mind, J? J+ j(s) may be rewntten in the form 

r 1 d 


w»> = /.’ ^ RHS]ds - 


(A.2) 


From (3.2a) we may find the dependence of v j+ ,(t + At) on the continuation parameter - 
(for simplicity we omit the explicit dependence on time): 


which in view of 
yields 


„ i+1 («) = \[v(s) + v j+1 ] - X[g i+ i - g{v{s))), 
J 2 £ 


-V(v i+i (a)) = + V(»W)1 ' A ”f+i- 

ds ^ 


In a similar manner, we have 


—U(v(s)) = -U'(v{s)) • Av i+ i, 
ds 


and Leibnitz rule gives us 

d 


At> i+ »- 


L\-\ r +l U'(v)g'(v)dv) - -At/>00)<7W) 
ds Jv{ 8 ) 

Substitution of (A.5), (A.6), and (A.7) into (A.2) yields 

H. + ,( 9 ) = -Av i+i [j\ + V(vM)l ' [WM) - P>, + j(.))l* 

Next, we use le continuation u(r,s) = rule) + (1 - *,♦. 1» (*•*») >" ° rd “ l ° 
the last difference on the right as 

fi d 


(A3) 

(AA) 

(A.5) 

(A.6) 

(A.7) 

(A.8) 


tf'(«W) - tf(», + j(«)) = f 0 


(A.9) 
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This equality comes about as follows: in view of (3.20b), (3.2a), v y+i (r,s) is given by 


W r ’ s ) = 2 ^ s) +t, ( r ’ 5 )] - %MM) - *(»(•))]; 

hence, v J+ ^(l,s) = v(s), v J+ i(0 ,s) = t/ y+ .(s) and (A.9) follows. 

Noting that 

d 

~ v {r,s) = - Av /+r s, 

then by carrying out the differentiaton on the RHS of (A.9) we obtain, 

^ U '(v jH (r,s)) = —U"(v i+ i(r,s)) . [I - A ff '(t,(r,s))] . s Av J+L . 


(A. 10) 


(A.ll) 


(A. 12) 


Substituting (A.9), (A.ll) and (A. 12) into (A.8), we will end up with the desired identity 
| o.ZZ 1 1 1 . 


We close this section with the proof of Lemma 3.8. 

The piecewise linear interpolant of the gridfunction, chosen in (3.24), 


has a fixed slope at each cell: 


g(V) = “ ”>) + » 


9'{v(r, s)) = $'(t/(s)) 


= Ag '+* 

A W 


(A. 13) 


(A. 14) 


From| A .H) and (3.22) we obtain that in the case of quadratic entropy function where 


*/♦*(») = j(A«y + j)< 




Av 


>+k. 


(A. 18) 


Moreover, the difference g(v) - /(„) between two neighboring values o, and v M , covers 
an area of size, 

> f v *+ 1 A 

A J v . ^ ~ K v )) dv = 2^+i + ftJAvy+j ~ A // +1 f(v)dv. (A. 16) 

Thus, in view of (A.15) and (A.16), the desired inequality, (3.26), boils down to 

[g j+ i + g;} • Av j+ i - A f f( v )dv + i(A ~~^) 2 - ^(Au i+ i) 2 < 0. (A.17) 

}+? 


-i 

2 


1 Ao. 

2<^ 

To verify the inequality (A.17), we recall that by (3.2a), (3.2b) we have 

1 

8A "* ' ' ' 2' m/ ' gA 


9rn + 2 )) + 8A ^ f(v m (t) - -f' m ) + m = jj + (A.18) 
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and Taylor’s expansion yields 


9m = /. + k(l - *P) + 0(A. J+i )*, p = MW‘»- 


(A.19) 


This enables us to write the first two terms on the left of (A.17) as 

fa ♦ •] • - * r Mdv - ^ ■ (Av,+i)2 + o(A ]f 0) 

Consider now the third term on the left of (A.17): by (A.19) we have 


i Av'. + i 

AA 9y+i = Wy+J + j(l - “0 ! ) • ' A * ,+ii 


(A.21o) 


inserting AA/ J+ r = + 0(An j+i )* into (A.21a), squaring the result and rearranging 


we obtain 


At/ i 

i(AA Sj+ .)’ = |(AV 4 ) 2 + |(1 - *f) ■ ' (A " i+4)!+ 


^ V j+k 


At;' 

L(1 _ 4/j*)* ■ (k 11 )’ • + °( Av >+i )S - 

28^ A «i+* 


(A.2lb) 
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We note that the cubic term on the right of (A.20), (A.21b), consists of the error in the 
trapezodial rule 

. /•«.•+. A 

A , 


j(/(v,+i) + /(®i)] A »i+l - A {,' f ^ dv ~ 12 / " t ” (l " ' 
as well as additional contributions which are of the same order of magnitude 

0(Av,- + 4 ) S < A • |/“(»W)1 • (Auy+t) 3 . 

Inserting (A.20), (A.21b), and (A.22) into tie inequality (A.17) gives us 


{A. 22) 


1-4/3 2 

8 


(A« J+ i) 5 


2Av J+ i Av ? +* 16 V Au i+ij 


+ 


+ A-max[/>(x))]-(A Vi+ ,) 3 <°- 

1 (>1.23) 

Finally, since ■} and ^ agree in sign with A.**, the expression inside the left brackets 

can be upper bounded by 


[•••]< 




l„ Av !>ii , Izjgli^li 




2' AVy + l 


16 Au ;+ i 


{AM) 
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1 ft hmitatl ° n ’ ^ < 2 > the Sum of the last two terms is nonpositive, and we 

left with the inequality 


are 


max 


v+i 


At v+* ’ 


- 1 


+ A • max[/>(x))] . Av j+k < o, 


which is met by the choice of entropy satisfying limiter in (3.25a), (3.25b). 
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